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We derive an analytic connection between the screened self-consistent effective potential from 
density functional theory (DFT) and atomic effective pseudopotentials (AEPs). The motivation 
to derive AEPs is to address structures with thousands to hundred thousand atoms, as given in 
most nanostructures. The use of AEPs allows to bypass a self-consistent procedure and to address 
eigenstates around a certain region of the spectrum (e.g., around the band gap). The bulk AEP 
construction requires two simple DFT calculations of slightly deformed elongated cells. The ensuing 
AEPs are given on a fine reciprocal space grid, including the small reciprocal vector components, 
are free of parameters, and involve no fitting procedure. We further show how to connect the 
AEPs of different bulk materials, which is necessary to obtain accurate band offsets. We derive a 
total of 20 AEPs for III-V, II- VI and group IV semiconductors and demonstrate their accuracy and 
transferability by comparison to DFT calculations of strained bulk structures, quantum wells with 
varying thickness, and semiconductor alloys. 
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I. INTRODUCTION 

The main motivation for deriving atomic effective 
pseudopotentials (AEPs) relics on the computational 
savings that may be gained by circumventing the self- 
consistent optimization of the density, which is at the 
heart of density functional theory (DFT). The aim is 
thereby to address the omnipresent nanoscale, where the 
relevant structures include thousands, and often hundred 
thousands, of atoms. The philosophy behind an AEP, 
not only bypasses the need for a self-consistent solution, 
but also allows to focus on a selected part of the eigen- 
value spectrum. The number of "bands" that scales with 
the number of atoms in DFT, is now independent of the 
system size. This latter characteristic is ideal for the 
study of optical properties or transport, where mainly 
the energy states around the band gap of the materials 
are involved. 

The idea to replace the strong Coulomb potential of 
the nucleus together with the bound electrons by an ef- 
fective, weaker (pseudo)potential is not newii 2 -. The first 
generation of the method^r— was based on the fitting of 
an empirical pseudopotential constructed to reproduce 
experimentally determined energy levels of bulk crystals. 
In this procedure, known as the empirical pseudopoten- 
tial method (EPM), the entire band structure as well as 
optical properties of bulk solids can be obtained with 
good accuracy with only very few Fourier components of 
the pseudopotential (in the order of three). These few 
parameters were simply used as fitting parameters and 
adjusted to reproduce known quantities. 

In the original EPM, the parameters were directly 
linked to the investigated bulk structure through the use 
of certain G-vectors and transferability (i.e. the applica- 
bility of one potential to different lattice constants or dif- 
ferent structures) was not supported by this type of con- 
struction. A simple extension via interpolation to neigh- 
boring points in G-spacO" 9 ^ or a full analytic dependence 



of the pseudopotential on the length of the G- vectors^— 
was a natural extension of the method. With these con- 
tinuous empirical pseudopotentials, together with mas- 
sively parallel codes^~— or with favorable basis setsi£~— , 
the electronic and optical properties of nanostructures 
could be accurately calculated^— , as well as transport 
propertie; 
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With the aim of getting an improvement in the quality 
of the wave functions, and hence of the transferability, a 
second generation of empirical pseudopotentials was de- 
veloped by Wang, Fu and Zunge r 28 ' 29 . In this approach, 
a series of bulk self-consistent DFT calculations, in the 
local density approximation (LDA), are performed using 
a few different structures (zincblcndc, wurtzitc, ...) at 
different lattice constants. The screened local effective 
potentials from DFT are transformed to G space and fit- 
ted by a set of Gaussian functions. This results in a con- 
tinuous representation of the atomic empirical potentials 
in G space with a good transferability between different 
bulk structures. However, as the structures used in the 
potential generation arc all bulk systems, no information 
on the long range (small G-vector) inter-atomic inter- 
action is available. Hence the potentials at G-vectors 
shorter than 2n divided by the lattice constant, remains 
unknown. Unfortunately, it is precisely this range of 
G-vectors that becomes relevant for the calculation of 
nanostructures, where the shortest relevant G-vectors is 
given by 2tt divided by the cluster size. In the original 
wor k 28 ' 29 , the Gaussian function was interpolated from 
the shortest known G-vector (2-k/o) to the G=0 point. 
The potential at the G=0 was adjusted to reproduce ex- 
perimentally known work functions^ and later using an 
atomic model potential and screening 2 ^. Furthermore, 
the energy cut-off was significantly reduced to save com- 
putational resources, so that an empirical Gaussian func- 
tion had to be introduced to compensate for this reduc- 
tion. Although the quality of the ensuing semiempirical 
pseudopotentials (SEPM) must be expected to be signif- 
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icantly higher than the EPM, these difficulties may have 
been the reason for the rather limited use of this potential 
generation. 

In this paper, we introduce a new generation of pseu- 
dopotentials following the philosophy of the EPM and 
SEPM to derive a pseudopotential that represent the 
scattering due to all the electrons. In this sense, it 
is distinct from the ionic pseudopotentials derived for 
DFT— ~— , although we use a separable formalism to treat 
the non-local part of the potential stemming from this 
community^. We derive our atomic effective pseudopo- 
tentials (AEPs) via an analytic connection to the effective 
crystal potential from DFT. The use of large and judi- 
ciously chosen slightly deformed supercells gives us access 
to the long range response of the potential. We therefore 
derive the pseudopotentials on a very dense grid of G 
vectors, including both, the long and short range inter- 
action regimes. The method is free of parameters and 
does not involve any fitting procedure and is as ab initio 
as the underlying DFT calculations. This simplicity al- 
lows us to generate, with a minimum of effort, unique (for 
a certain DFT norm-conserving pseudopotential) AEPs. 
We have successfully generated AEPs for 20 compound 
semiconductor materials. A comparison between the re- 
sults of our code LATEPP^ using our AEPs and DFTi2. 
shows differences for the bulk calculations in the order 
of a few tenth of meV, which originates from our spheri- 
cal approximation of the local part of the potential. We 
have tested the transferability of the potential on dif- 
ferent types of quantum wells and find very satisfying 
results, demonstrating a high level of transferability. 

II. METHODOLOGY 

Our methodology is based on DF T 41 i 42 where the cen- 
tral quantity is the effective Kohn-Sham potential V4ff 
that describes all the interactions of a single electron with 
its environment: 

V oS (r) = V cxt (r) + ^ Hartro Xr)] + V xc [n(r)} (1) 

occ 

n ( r ) = ^2\^( r )\ 2 > 

i 

i 



This way, the ionic pseudopotential has been separated 
into a local /-independent part and non-local operator. 



where n(r) is the charge density of all occupied single- 
particle states ipi . The second key component is given by 
the transformation from the true Coulomb potential 



7 p 2 

V^(r) ee V^(r; {R a }) = £ (2) 



into a pseudopotential V psp . The core electrons are 
thereby frozen in their atomic configuration following 
the frozen-core approximation^. The pseudopotential 
is usually written as a sum of angular momentum pro- 
jectors: 



impedes N a 

y P s P = J- £^^(|r-r Qn |)P ;m , (3) 

a n Ira 



where a describes the atom type and runs from one to 
-Species and n describes the atom number running from 
one to the number of atoms of type a, N a . The atomic 
positions are given by T an and P; m is the projection oper- 
ator. Here we note that all the pseudopotentials describe 
correctly the long-range Coulomb interaction, beyond a 
certain cut-off: 



«oi(r) = -— for r>rSf , (4) 



where Z a is the charge of the valence of the respective 
atom. It is therefor customary to perform the following 
transformation: 



(5) 
(6) 



The local potential must decay according to Eq. but 



^p Sp = £^ps Pl ioc ( | r _ Tan]) p lm + J2J2Ki(\r - r an \) - v™> loc (\r - r an \)}P t 

an lm an lm 

= E^ SP ' l0C d r - T »l)E^ +EEKl(k - Tanl) " < SP ' l0C (|r ~ T Q „|)]P 
an lm an lm 

= E< p ' loc (|r-T Q „|)+EE,5 UQ; (|r-T ml |)P iln . 
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is otherwise arbitrary. A good choice is to take one of 
the /-channels of the atomic pseudopotential. Since the 
local potential and all the /-components of the potential 
have the same long-range behavior, this formal transfor- 
mation leads to /-dependent potentials that differ from 
zero only within spheres with radius . This is numer- 
ically a decisive advantage. Atomic pscudopotcntials are 
constructed up to a certain maximum angular momen- 
tum Z m ax- From Eq. ((6]) it is clear that the truncation to 
( max only applies to the non-local parts of the potential. 
One could therefore argue that higher angular momen- 
tum component are incorporated in the local potential. 
This fact is not effectively used in the construction of 
ionic pseudopotentials but will have some importance for 
the AEPs we will derive. 

For the implementation of the projectors we follow the 
separable form given by Kleinman and Bylander— : 

V KB =V^ + Y,\xL B )Er{xL B \ , (7) 

l m 

with the Klcinman-Bylander eigenvalue given by: 



e kb = (SViui\SViui) 
1 {u t \5Vi\ui) 



(8) 



and Xhn t ne normalized Klcinman-Bylander projectors 

\6Vl4>l m ) 



xL u ) 



((mSVilSVim)) 1 /* 



(9) 



The effective crystal potential from Eq. ([T]) is now given 
more specifically as: 

yeff = yp8 p,loc + F Hartree[ n ] + F xc[ n j + ^ 

J2\xL B )Er(xf n B \ ■ (11) 
in i 

During the self-consistent cycle of the Kohn-Sham equa- 
tions (Eq. ([I])) the density is updated until the ground 
state density n is found. The starting point for the 
derivation of our AEPs is the local part of the self- 
consistent effective potential: 



F loc,eff (r) = F p SP ,loc (r) + ^^cfj + ^ 



scf "1 



with 



V psp ' loc (r)=J2< P ' loC (\r-r an \) 



(12) 



(13) 



where tiP s P- loc is the local potential we defined in Eq. ([5]) 
and is an output of the pseudopotential construction 
procedure (e.g. after Hamann, Schliiter and Chiang^, 
Kerker— , Vanderbilt^I or Troullier and Martini!) and 
usually includes one of the angular momentum com- 
ponent of the ionic pseudopotential. The Hartree and 
exchange-correlation parts are, however, non-spherically 
symmetric. 



From the self-consistent real space potential we can 
move to reciprocal space using: 
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We rewrite the local effective potential as a sum of atom- 
ccntcrcd potentials: 

Nspccics N a 

y iocc ff(r)= J- 5> Q (r-r Q „) . (15) 

a n 

With this transformation we can rewrite the total re- 
ciprocal space potential (via variable substitution) as a 
Fourier sum: 

^Species N a 



V 



l0 c,c ff(G) = _i_ ^ y: 



-iG-T a n 



v a {G) (16) 



a n 



with 



v a (G) = I v a (r)e~ lGr d 3 r . (17) 



III. AEP CONSTRUCTION PROCEDURE 

A. Binary systems 

The AEPs of a binary system cannot be directly ex- 
tracted from the screened local effective potentials of 
DFT, but can be determined through v + and w_, 

V+ = Va+V c , (18) 
V- = V a -V c , 

where v a {v c ) stands for the anion (cation) potential. In 
the following we will use a binary system as prototype, 
but the procedure can be simplified to pure material by 
setting v a = v c . We start by explicitly dividing Eq. ([15]) 
into cation and anion components: 



N,. 



^c(r) = E^( r -^)+E^( r -^) > ( 19 ) 

t=l 3=1 

with N a (N c ) as the number of anion (cations) in the 
structure. 

Keeping the same atomic structure, and merely inter- 
changing the atom types (index i and j), we rewrite: 



N c 



V^(r)=Y,v c (r-7i) + ^2v a (r-^) . (20) 
i=i i=i 

By adding and subtracting V^J and , we find the 
following expressions for v + and V- : 

N atoms 



V& +2 \r) = vU(r) + vM(r) = £t, + (r-r n ) , (21) 

N atom s 

V&- 2) (r) = vg>(r)-Vg>(r) = £ (-l)' 1+ V(r - r„|22) 
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The Fourier transform of the last two equations leads to 
the relations: 



V, 



(1+2) 



loc 



(G) 



V, 



(1-2) 



loc 
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i 
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e iG-rv(l+2) 



E 

n=l 

AT otom , 

E 



(r)dr 
v+(G) 



JG-T„ 



(-1) 



(23) 

V-(G) ( 24 ) 



The potentials v + (G) and V—(G) are complex, but into 
the scheme of the spherical approximation used in this 
work, only the real part of the atomic potentials is of 
interest. We can extract the spherically averaged AEPs 
from Eqs. (|2"3"]) and ([2~4")l . by separating them into their 
real and imaginary components, and write: 



4 A 



n 



(\G\)=Re[v ± (G)] = 
Re[v£ c ±2) {G)] , Im[vi\ ±2 \G)] 



" loc 

Pi 



P±a± 
HI + a\ 



(25) 
(26) 



where 



N ato 



8+ = sin ( G ' T ' 1 ) > 

n=l 

N atom s 

a+ = ^ C0S ( G ' T n) , 

n=l 

Natoms 

P-= E (-l)" +1 sin(G-r„) 

n=l 

a_= X] (-l) n+1 cos(G-r„) 



(27) 



(28) 
(29) 



n=l 



With 4 A (G), we can use Eqs. (fT8|) to obtain the 
atomic potentials w a and v c in G-space. 

The previous development can be used with any bi- 
nary system. However, not every structure is suitable to 
extract the AEPs. It is, for example, of interest to obtain 
the potential for a high enough density of points, not only 
because that leads to the inclusion of the long range in- 
teraction, but also because that allows the utilization of 
a simple interpolation of points as a mean of constructing 
a continuous representation of the potentials in G-space. 

To generate the atomic potentials, we used a 
zincblende structure made of 24 atoms and elongated 
along the [100] direction, regardless of the ground state 
structure of the different binary systems (wurtzite (wz) , 
zincblende (zb) or rocksalt (rs)). We apply a compressive 
and a tensile deformation of 5% along the slab direction 
in order to break the symmetry of the crystal and allow 
the extraction of the long range interaction. Note that 
a large supercell with ideal bulk atomic positions leads 
to a fine G-vector mesh, but the structure factor is ex- 
actly zero at all the points, except at the G- vectors of 



the bulk unit cell. As it should be, the large periodic 
supercell carries the same information as the unit cell, 
hence the need for the deformation. 

In Eqs. (|2"3"]l and ([2~4")l we have assumed that the atomic 
potentials are independent of the atomic positions, which 
is strictly true for the bulk unit cell calculation but only 
approximately for the slab calculation. The main pur- 
pose to use the deformed slab calculation is to obtain 
points for small G-vectors that are unavailable from a 
bulk unit-cell calculation. The slab calculations deliver, 
however, results for the entire range of G-vectors and 
these results are identical (up to numerical accuracy) 
with the results obtained from a bulk unit-cell calcula- 
tion for the large G-vectors. For the shortest available 
G-vector from a unit-cell calculation, G c , we observe a 
small deviation between the slab and the uni-cell results. 
We therefore correct the slab results at this specific G- 
vectors by a Gaussian correction to exactly obtain the 
unit-cell result (where the assumption of unique atomic 
potentials is exactly fulfilled). The shape of the Gaussian 
correction is fully defined by assuming that the correc- 
tion has dropped to 1% of its original value at G equal 
to two times G c . This procedure is not strictly necessary 
but slightly improves the results and requires only little 
effort. 

The orientation of the slab is chosen along the [100] 
direction because it allow for the extraction of w^ A (|G|). 
This becomes clear if we write the local potential as 



N atom3 /2 

Vi oc (r) = ^2 {Vair-Rn + r^ + Vcir-Rn-rn)} , 

n 

(30) 

with r n the position of the atom n measured from the 
anion-cation bond center R n . In reciprocal space we ob- 
tain, 

N atoms /2 

Moc(G) = 75 E e lGRn {e~* Grn Va(G) + e* Gr -v c m} 



Natoms/2 



1 J- e* GR "{cos(G-r n )v + (G) 



(32) 



;sin(G-r r >-(G)} 



We see, in the last equation, that whenever G and r 
are perpendicular (as can happen in the [110] orientation, 
where the cell can be written by placing the atoms in 
planes perpendicular to the direction of the structure) 
the quantity multiplying v—(G) becomes zero and u_(G) 
disappears from the equations. On the other hand, this 
expression shows us that V- (G) becomes less and less 
important as G approaches zero. 

To summarize, the generation of AEPs comprises four 
steps: i) The DFT-LDA calculations of two 24 atoms 
distorted cell structures, where the atomic positions are 
interchanged, ii) The extraction of the effective atomic 
potentials from the screened local potentials obtained in 
step (i), by applying Eq. (|2"6")l for G-vectors along the 
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slab direction leading to a high density of G-points. hi) 
Adjustment to the bulk properties in such a way that 
v± A (|G|) passes exactly trough i>:| A ( | G c | ) . In this step 
a DFT calculation of the bulk system is needed, iv) A 
splines interpolation of the calculated points, in order to 
have a continuous form of the potentials in G space. 



B. Heterostructures and band-offsets 

Until now, the AEPs were derived from DFT calcula- 
tions of periodically repeated simulation cells of one ma- 
terial. For periodic systems, the average effective crystal 
potential, which corresponds to the V(G = 0) compo- 
nent of the potential is ill defined and fixed to an arbi- 
trary values in DFT codes. This means, that the ener- 
getic position of the eigenstates can be arbitrarily shifted 
along the energy axis. This also means that the eigen- 
values of two DFT calculations for different materials, 
or even different structures of the same material, can- 
not be compared on an absolute scale. Only differences 
of eigenvalues within one DFT calculation can be for- 
mally compared (see, e.g,. RefliU for a discussion). In a 
heterostructure made of two materials A/B and treated 
with periodic boundary conditions, the V(G = 0) com- 
ponent remains ill defined and is arbitrary. However, the 
difference between the "local" average potential in the 
A and in the B regions, as far as this quantity can be 
meaningfully defined and calculated, represent the band- 



offset and is a very real and important quantity. The 
sign of the band-offset between materials will lead the 
the localization of the carriers in either the A- or B- 
part of the structure. Calculations of band offsets based 
on DFT— ~ involve the calculations of heterostructures 
and is a rather subtle task. 

We incorporate the effects of band-offsets in our AEPs 
by linking them together using DFT calculations of het- 
erostructures. The procedure involves the interchange of 
the cation anion positions within a QW heterostructure 
formed by two slabs of different materials (A and B). The 
local potentials can be written as 



N a 








-n)- 




t=i 












£<0 


- Tfe) 




k=l 




i=i 


N? 




N a 


£0 


~Ti)- 




i=l 




j=l 








k=l 


-r k ) 


1=1 



Extracting the real part of the potential for slab B (v±) 
we obtain two equivalent expressions 



Re\v?(G) 



^Re[v[£ 2 \G)]-^Re[v£(G)] 



n 



Re[v«(G)] = j^Im[v[lf 2) 
P± 



(G)} 



0± 



%Re[u±(G)] 



^I m [vi(G)} + %im[vi(G)} 



B 



0! 



Im[v±\ 



G)] 



Im[v£(G)] 



(35) 
(36) 



where (3± and a± are given by Eqs. (|29|) . but restricted 
to each one of the slabs A oi B. 

We therefore have a set of two equations and four un- 
knowns. Fortunately, the imaginary components, which 
are basically the result of the deformation or loss of 
sphericity of the potentials due to the environment, van- 
ish in most of the points and the expressions simplify 
to 

Re[Vi{G)] = —Re[v£ c ±2) (G)] - %i?e[V^(G)] (37) 
ct_l_ ct_l_ 

Re[V£{G)\ = ^Im[vZ ±2 \G)] - ^Re[v£(G)] (38) 

P± P± 

We now fix the AEPs of slab A to the potential derived 
from the bulk calculation and extract from either one of 
the equations the potentials for slab B. 



In this paper, we link the potentials by fixing each po- 
tential once it has been linked and subsequently using it 
as a basis to connect the new potential (linking B to A, 

C to B, D to C, A <- B <r- D < ). We defined 

the linking order by using QWs with the lowest possible 
lattice mismatch among the studied materials. We used 
a 48 atoms zincblende QW structure, along the [100] di- 
rection, made of two slabs (A and B) of equal sizes, with 
a lattice constant equal to the average of the lattice con- 
stants of the two materials, including the interfaces, and 
without any relaxation. 
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TABLE I: Absolute volume deformation potential (in eV) for 
the valence band maximum (ajf sjv/ ). 



T T T \ J 

ill- V 


VBAI 

a v 


T T ^ I T 

11- VI 


VBAI 

a v 


1V-1V 


VBAI 

a v 


InSb 


1.61 


CdTe 


0.89 


Ge 


2.23 


AlSb 


2.65 


CdSe 


0.90 


SiGe 


2.30 


InAs 


i 7n 


Znle 


0.99 


Ol 


2.38 


GaSb 


1.91 


ZnSe 


1.23 


SiC 


4.30 


InP 


1.83 


MgO 


1.66 


c 


2.16 


AlAs 


2.93 


ZnO 


0.48 






GaAs 


2.24 










A1P 


2.60 










GaP 


1.98 











C. Deformation Potentials 

The v+{G = 0) component of the potential does only 
shift the eigcnspectrum along the energy axis and has no 
meaning as it comes out of a DFT calculation. However, 
this value influences the conduction and valence band 
deformation potentials individually. Again, the gap de- 
formation potential is (as it should) independent of this 
value. We decide to fix the v + (G = 0) to reproduce the 
DFT values of the valence band deformation potential 
according to Li et aL— ^ for most of the materials. For 
SiGe, we averaged the values of pure Si and pure Ge. 
For SiC, we reproduced the absolute deformation poten- 
tial reported by Lambrecht et al^. In the case of MgO, 
we used the value reported by Zhu et aL— The values we 
used are given in Table HI 



IV. ASSESSMENT OF THE POTENTIAL 
QUALITY 

A. Used norm-conserving pseudopotentials 

We started by generating a set of norm-conserving 
pseudopotentials to be used in the DFT calculations. 
The pseudopotentials used to calculate the screened lo- 
cal potentials in DFT are directly related to the gener- 
ated AEPs. The AEPs have to be used together with 
the non-local part of the norm-conserving pseudopoten- 
tials in any future calculation. As our main interest in 
extracting AEPs is the reduction in computational time 
when studying semiconductor nanostructures, we used 
soft pseudopotentials generated with the FHI code^, by 
including only up to p electrons when possible, and go- 
ing up to d electrons when strictly required. In the case 
of the group III-V and group IV semiconductors, we al- 
ways used the default FHI values for the real space cut-off 
radius of the different angular momentum channels. In 
the case of the II- VI compounds, we increased the real 
space cut-off radius as it leads to faster convergence with 
the enegy cut-off at a very modest cost in accuracy (see 
later) . 



B. Comparison with DFT for bulk: errors induced 
by the spherical approximation 

The first test of our potentials is against the bulk re- 
sults obtained by a self-consistent DFT calculation nat- 
urally performed without spherical approximation. We 
refer to these calculations as "DFT" in the first column 
in Table |TT] where we report the value of the band gaps. 
In the second column, "AEP" , we report the results us- 
ing our AEPs. The third column gives the energy dif- 
ference between both approaches and is the result of the 
errors stemming from the spherical approximation. The 
last column shows the DFT relaxed lattice constants (a ) 
used to generate the atomic potentials. We have, in gen- 
eral, errors in the energy band gap of few tenths of meV. 
With the exception of carbon, with an error of 160 meV, 
which is still small compared to the LDA band gap of over 
5 eV. In Table ILT1 we find the first tests of the transferabil- 
ity of some AEPs. In the case of CdSe, MgO and ZnO, 
we have calculated the bang gap of stable structures dif- 
ferent than zincblende. We find that the deviations due 
to the spherical approximation is in the same order of 
magnitude as the deviations we obtain by changing the 
crystal structure. Note that InAs and InSb have a van- 
ishing band gap in LDA, so that we chose to give the gap 
at the L-point in Table |Hj 



C. External and Internal Strain Situations 

In the following test, we selected one III-V, one II- VI 
and one group IV material, namely A1P, ZnSe and SiC, 
and studied the band gaps under different strain condi- 
tions. Wc start by studying the bulk materials under the 
effects of uniform expansion and compression and show 
the results in Fig. [TJ We see a quasilinear dependence 
in all cases, going from large band gaps, when the sam- 
ple is compressed, to smaller band gaps in the expanded 
situation. The AEP results follow well the self consis- 
tent DFT calculations. In Fig. [2] we plot the deviations 
between the AEP and the DFT results. This deviation 
becomes important when the effective potential starts to 
be significantly modified by self-consistent effects. Until 
around 5% the approximation of rigid overlapping poten- 
tials remains accurate. The value of 5% is very large for 
a hydrostatic deformation and will not be encountered in 
any experimental situation. 

Variations of the lattice constant are translated in 
terms of variations in the reciprocal space mesh. If the 
lattice constant is expanded, the G-space mesh becomes 
denser and the AEP calculation probes G-vectors that 
are shorter than the critical length G c defined by the 
relaxed lattice. The calculations done with compressed 
lattice constants accordingly probe G-vectors longer than 
G c . This constitutes a good indication of the transfer- 
ability of our AEPs. 

Fig. [5] also shows the deviations of the band gap ener- 
gies when the calculations are done with the AEPs before 
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TABLE II: Comparison of the single-particle band gap gap 
at the T-point or at the L-point when indicated by (*). The 
bulk crystals have mostly the zincblende structure and, in a 
few cases, the rocksalt (rs) or wurtzite (wz) structures. For 
the wurtzite structure, both lattice parameters, c and a, are 
given. 



Compound 


JJr 1 (ev ) 
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FIG. 1: Bandgap for selected bulk materials, as a func- 
tion of the lattice constant deformation calculated using self- 
consistent DFT and using our AEPs. 
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FIG. 2: Deviations of the AEP calculations from the DFT 
results for different bulk materials as a function of lattice de- 
formation. 



and after ("link" notation in Fig. [2]) the linking procedure 
(Sec. IIIIBj) . As expected, we see that the results of the 
AEP calculations are identical in the compressed situa- 
tions. We have differences between the results of only 
few meV in the expanded condition, which is acceptable. 

One more test performed in this section is aimed at 
testing not only the response of the material to strain, 
but also the long-range quality of our potential. For this 
purpose we constructed a 32 atom and a 36 atom cell 
elongated along the [110] and the [100] directions, respec- 
tively We compress one half of the cell and expanded the 
other half, along their respective orientations. The lat- 
tice constant in the perpendicular direction was kept at 
the bulk value. With this construction, the interface be- 
tween the expanded and the compressed regions has the 
ideal lattice constant as well, and the length of the cell 



remains constant for different deformations. 

In Tables IIIII and IIVI we report the difference in the 
band gap between the AEP and the DFT calculations. 
As these slabs are binary systems, we can use again the 
potentials as originally derived from the deformed struc- 
tures or the atomic potentials obtained after the linking 
procedure. In all cases the variations of the band gaps 
with the deformations are very well reproduced by the 
AEP calculations. However, with a slightly better agree- 
ment for the [110] slab. The deviations of the binary 
potentials of AIP and ZnSe (second and third columns in 
Table IIII[) change by only few meV as the deformation is 
increased. The better quality in the results of the [110] 
sample may be attributed to the fact that the potential 
V- does not play any role in this structure. In the case of 
SiC, the results arc deteriorated by the linking procedure, 
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TABLE III: Energy band gap deviations of the AEP calcula- 
tions from the DFT results AE 9 = E 9 (AEP) — E 9 (DFT). The 
AEPs are obtained from the binary system (Binary) and after 
the linking procedure (Linked). The deformation is applied 
along the [110] direction. 



Deformation 


AE 9 (meV) 


A1P 


ZnSe 


SiC 


Binary 


Linked 


Binary 


Linked 


Binary 


Linked 


0% (Bulk) 


55.2 


55.2 


-18.6 


-18.6 


63.1 


63.1 


1% 


84.4 


82.3 


-20.7 


-25.2 


11.4 


-34.1 


2% 


81.6 


76.1 


-27.0 


-39.1 


-12.1 


-118 


3% 


81.4 


72.0 


-32.7 


-52.1 


-33.5 


-208 


4% 


83.4 


69.5 


-36.2 


-62.1 


-52.7 


-304 


5% 


87.8 


68.7 


-36.8 


-68.7 


-69.8 


-405 


6% 


94.6 


69.7 


-34.9 


-72.2 


-84.6 


-511 



TABLE IV: Equivalent to Table iHll for a deformation applied 
along the [100] direction. 



Deformation 


AE 9 (meV) 


A1P 


ZnSe 


SiC 


Binary 


Linked 


Binary 


Linked 


Binary 


Linked 


0% (Bulk) 


55.2 


55.2 


-18.6 


-18.6 


63.1 


63.1 


1% 


81.7 


72.5 


-22.0 


-26.5 


8.30 


-76.9 


2% 


67.7 


45.1 


-32.6 


-45.4 


-36.1 


-237 


3% 


52.6 


16.5 


-44.1 


-64.9 


-82.4 


-398 


4% 


37.8 


-11.4 


-54.9 


-82.6 


-128 


-558 


5% 


23.6 


-38.8 


-63.7 


-97.3 


-173 


-717 


6% 


9.02 


-66.5 


-71.7 


-110 


-217 


-793 



which is due to the large lattice mismatch, of around 20%, 
with respect to Si (Table [TTJ) . We therefore have to use 
a strongly deformed structure to connect the potentials 
of SiC to Si and C, which compromises the quality of 
the results. In this case, the quality of the linked AEPs 
could be improved by selecting a more suitable structure 
to link the potentials. 

Since the deviations between the DFT and the AEP 
results are only slightly different for the two different ori- 
entations, we conclude that the long-range response of 
the potential to the strained situation is not significantly 
affected by the nature of the interface. In the case of 
a [110] structure, the interface layers are composed of 
mixed anion-cation atoms, while the interfaces in a [100] 
structure are composed of pure anion and pure cation 
atoms; qualitatively very different situations. This rep- 
resents an important conclusion of our work. 



D. Semiconductor alloys: errors induced by the 
local average 

Here we considered 64 atom supercells of Ga^Ali-^As 
and GaaJni-zP alloys in the zincblende structure. The 
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FIG. 3: DFT and AEP band gap calculations for a 64 atom 
alloy system at different concentrations. 

alloying is simply done as a local average, where the AEP 
for an As atom surrounded by n Ga and (4 — n) Al atoms, 
is written as a weighted sum of the type: 

v As {nGa, (4 - n)Al} = ± („ v^ As + (4 - n) vf As ) , 

where v A iAs(GaAs) represents the AEP of As derived from 
the AlAs* (GaAs) structure. Note that the AEPs vf- As 
and v^ As should be different since the As AEP con- 
tains informations on the bonding and it's environment. 
In practice both quantities are surprisingly similar. In 
Fig. [3] we show the DFT results along with the results of 
the AEPs in a similar way as done before. The atomic 
concentration has been modified in such a way that we 
have a successive transition from one bulk material to 
the other (x = or x ~ 1). The atom positions have 
been relaxed in order to approach a realistic structure, 
by using the valence force field metho d 20 i 55 . 

The AEP and DFT calculations show the same quali- 
tative behavior and differ quantitatively by a magnitude 
similar throughout the composition range. These results 
demonstrate the validity of the linking procedure. It is 
worth mentioning here that AEP calculations using un- 
linked potentials lead to very different results as the band 
offsets are arbitrary. 

E. Quantum wells: errors induced by a lack of 
charge transfer at the interface 

The next test we perform is aimed at testing the qual- 
ity of the long-range response of the AEPs due to the 
presence of a heterostructure interface. In Fig. @] we 
plot the band gap as a function of the well (or barrier) 
thickness for a 36 atom [100] oriented A/B superlatticc 
with zincblende structure. As test systems we select: 
GaP/AlAs, CdSe/ZnTe and Si/Ge. The structures have 
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FIG. 4: DFT and AEP calculations for superlattices with 36 
atoms along the [100] direction. The size of the wells and 
barriers are modified trough the variation of the proportion 
of the compound materials. 



been relaxed by using and averaged lattice constant per- 
pendicular to the superlattice direction and according to 
Poisson's ratio along the superlattice. 

Through this configuration procedure, the structures 
all have different dimensions, this means that they re- 
quire a different set of G- vectors in reciprocal space. Be- 
sides, the G mesh is rather dense, due to the size of the 
sample. This represents a highly demanding test for our 
AEPs since the full range of G-space is required. Fur- 
thermore, in each one of these calculations a consider- 
able number of potentials is needed due to the presence 
of the interface. For instance, for the GaP/ AlAs super- 
lattice, we need eight AEPs; four for each one of the slabs 
(GaP and AlAs) and four for the two interfaces (A1P and 
GaAs). 

All the AEP calculations show good agreement with 
the DFT calculations in Fig. H The CdSe/ZnTe super- 
lattice shows less deviations in the energy band gap, but 
the results are also less parallel to the DFT calculations 
in comparison to the GaP/AlAs and Si/Ge cases. This is 
due to the increase in the real space cut-off radii applied 
to the norm-conserving pseudopotentials of the II- VI ma- 
terials, which leads to a loss in transferability in both, 
DFT and AEP calculations. However, the deviations for 
the CdSe/ZnTe superlattice never exceed 30 meV, well 
within our target accuracy. 

In order to see the effects on the quality of the wave 
functions, we selected the GaP/ AlAs supercell with 14 
layers of GaP and 4 layers of A1P and show in Fig. [5] 
the square of the VBM and CBM wave functions after 
averaging over the (100) planes. The figure shows a type 
II superlattice, where the VBM wave function is mostly 
localized in the GaP region and the CBM wave functions 
in the AlAs region. We compare the AEP with the DFT 
state densities and notice a good agreement for the enve- 
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FIG. 5: Square of the CBM and VBM wave functions in a 
[100] GaP/ AlAs superlattice with 14 layers of GaP and 4 
layers of A1P. The wave functions have been integrated in 
planes along the axes of the superlattice. The calculations 
are done using DFT and our AEPs. 



lope of these strongly oscillating functions, but also some 
differences in the values of maxima and minima within 
these atomic oscillations. Looking at the CBM densities 
carefully, we notice a slight increase of the density in the 
AlAs layer compared to the AEP case. A slightly less 
confined CBM wave function is in agreement with the 
small energetic differences found in Fig. 0] for this case. 
Of relevance for the test is the fact that the long range 
response of the potentials are called upon to describe the 
long-range oscillations of the envelope functions and the 
asymmetry on the localization of the wave functions due 
to the different character of the interfaces, GaAs on one 
side and A1P on the other. 



V. SUMMARY 

In summary, we have developed a methodology to de- 
rived atomic effective pseudopotentials (AEPs) from sim- 
ple DFT calculations. The procedure involves DFT cal- 
culations for two slabs with compressed and expanded 
regions. We establish an analytic connection between the 
DFT effective Kohn Sham potentials and the AEPs. The 
procedure is therefore free of parameters and does not 
involve any fitting procedure, which represents the main 
achievement of this work. Wc furthermore establish a 
method to connect the AEPs for different materials, rig- 
orously from DFT calculations. Our AEPs intrinsically 
contain band-offsets between materials that are inherited 
from the DFT calculations, without having to explicitly 
calculate them. We demonstrate the accuracy and trans- 
ferability of the AEPs for a total of 9 III-V, 6 II- VI and 
5 group IV elements and find very good agreement with 
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self-consistent DFT calculations. 
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